{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variance Ratio Test - Lo and MacKinlay (1988)\n",
    "\n",
    "refer to https://mingze-gao.com/measures/lomackinlay1988/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LoMacKinlay.py\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import pandas as pd\n",
    "from pprint import pprint\n",
    "\n",
    "name = 'LoMacKinlay1988'\n",
    "description = 'Variance ratio and test statistics as in Lo and MacKinlay (1988)'\n",
    "vars_needed = ['Price']\n",
    "\n",
    "\n",
    "@jit(nopython=True, nogil=True, cache=True)\n",
    "def _estimate(log_prices, k, const_arr):\n",
    "    # Log returns = [x2, x3, x4, ..., xT], where x(i)=ln[p(i)/p(i-1)]\n",
    "    rets = np.diff(log_prices)\n",
    "    # T is the length of return series\n",
    "    T = len(rets)\n",
    "    # mu is the mean log return\n",
    "    mu = np.mean(rets)\n",
    "    # sqr_demeaned_x is the array of squared demeaned log returns\n",
    "    sqr_demeaned_x = np.square(rets - mu)\n",
    "    # Var(1)\n",
    "    # Didn't use np.var(rets, ddof=1) because\n",
    "    # sqr_demeaned_x is calculated already and will be used many times.\n",
    "    var_1 = np.sum(sqr_demeaned_x) / (T-1)\n",
    "    # Var(k)\n",
    "    # Variance of log returns where x(i) = ln[p(i)/p(i-k)]\n",
    "    # Before np.roll() - array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])\n",
    "    # After np.roll(,shift=2) - array([8, 9, 0, 1, 2, 3, 4, 5, 6, 7])\n",
    "    # Discard the first k elements.\n",
    "    rets_k = (log_prices - np.roll(log_prices, k))[k:]\n",
    "    m = k * (T - k + 1) * (1 - k / T)\n",
    "    var_k = 1/m * np.sum(np.square(rets_k - k * mu))\n",
    "\n",
    "    # Variance Ratio\n",
    "    vr = var_k / var_1\n",
    "\n",
    "    # a_arr is an array of { (2*(k-j)/k)^2 } for j=1,2,...,k-1, fixed for a given k:\n",
    "    #   When k=5, a_arr = array([2.56, 1.44, 0.64, 0.16]).\n",
    "    #   When k=8, a_arr = array([3.0625, 2.25, 1.5625, 1., 0.5625, 0.25, 0.0625])\n",
    "    # Without JIT it's defined as:\n",
    "    #   a_arr = np.square(np.arange(k-1, 0, step=-1, dtype=np.int) * 2 / k)\n",
    "    # But np.array creation is not allowed in nopython mode.\n",
    "    # So const_arr=np.arange(k-1, 0, step=-1, dtype=np.int) is created outside.\n",
    "    a_arr = np.square(const_arr * 2 / k)\n",
    "\n",
    "    # b_arr is part of the delta_arr.\n",
    "    b_arr = np.empty(k-1, dtype=np.float64)\n",
    "    for j in range(1, k):\n",
    "        b_arr[j-1] = np.sum((sqr_demeaned_x *\n",
    "                             np.roll(sqr_demeaned_x, j))[j+1:])\n",
    "\n",
    "    delta_arr = T * b_arr / np.square(np.sum(sqr_demeaned_x))\n",
    "\n",
    "    # Both arrarys are of length (k-1)\n",
    "    assert len(delta_arr) == len(a_arr) == k-1\n",
    "\n",
    "    phi1 = 2 * (2*k - 1) * (k-1) / (3*k*T)\n",
    "    phi2 = np.sum(a_arr * delta_arr)\n",
    "\n",
    "    # VR test statistics under two assumptions\n",
    "    vr_stat_homoscedasticity = (vr - 1) / np.sqrt(phi1)\n",
    "    vr_stat_heteroscedasticity = (vr - 1) / np.sqrt(phi2)\n",
    "\n",
    "    return vr, vr_stat_homoscedasticity, vr_stat_heteroscedasticity\n",
    "\n",
    "\n",
    "def estimate(data):\n",
    "    \"A fast estimation of Variance Ratio test statistics as in Lo and MacKinlay (1988)\"\n",
    "    # Prices array = [p1, p2, p3, p4, ..., pT]\n",
    "    prices = data['Price'].to_numpy(dtype=np.float64)\n",
    "    result = []\n",
    "    # Estimate many lags.\n",
    "    for k in [2, 4, 6, 8, 10, 15, 20, 30, 40, 50, 100, 200, 500, 1000]:\n",
    "        # Compute a constant array as np.array creation is not allowed in nopython mode.\n",
    "        const_arr = np.arange(k-1, 0, step=-1, dtype=np.int)\n",
    "        vr, stat1, stat2 = _estimate(np.log(prices), k, const_arr)\n",
    "        result.append({\n",
    "            f'Variance Ratio (k={k})': vr,\n",
    "            f'Variance Ratio Test Statistic (k={k}) Homoscedasticity Assumption': stat1,\n",
    "            f'Variance Ratio Test Statistic (k={k}) Heteroscedasticity Assumption': stat2\n",
    "        })\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'Variance Ratio (k=2)': 1.0003293867428107,\n",
      "  'Variance Ratio Test Statistic (k=2) Heteroscedasticity Assumption': 0.00032904650491573964,\n",
      "  'Variance Ratio Test Statistic (k=2) Homoscedasticity Assumption': 0.3293865781172764},\n",
      " {'Variance Ratio (k=4)': 1.0007984480057008,\n",
      "  'Variance Ratio Test Statistic (k=4) Heteroscedasticity Assumption': 0.000425953554365409,\n",
      "  'Variance Ratio Test Statistic (k=4) Homoscedasticity Assumption': 0.4267881978179488},\n",
      " {'Variance Ratio (k=6)': 0.9999130202975436,\n",
      "  'Variance Ratio Test Statistic (k=6) Heteroscedasticity Assumption': -3.511758587335342e-05,\n",
      "  'Variance Ratio Test Statistic (k=6) Homoscedasticity Assumption': -0.03518500446740915},\n",
      " {'Variance Ratio (k=8)': 1.0001094011344323,\n",
      "  'Variance Ratio Test Statistic (k=8) Heteroscedasticity Assumption': 3.69227065980853e-05,\n",
      "  'Variance Ratio Test Statistic (k=8) Homoscedasticity Assumption': 0.03698431520284624},\n",
      " {'Variance Ratio (k=10)': 1.0007024101299271,\n",
      "  'Variance Ratio Test Statistic (k=10) Heteroscedasticity Assumption': 0.00020772753506398227,\n",
      "  'Variance Ratio Test Statistic (k=10) Homoscedasticity Assumption': 0.20803582207648225},\n",
      " {'Variance Ratio (k=15)': 1.0022173139633859,\n",
      "  'Variance Ratio Test Statistic (k=15) Heteroscedasticity Assumption': 0.000521307044544808,\n",
      "  'Variance Ratio Test Statistic (k=15) Homoscedasticity Assumption': 0.5219816274022102},\n",
      " {'Variance Ratio (k=20)': 1.003804866170505,\n",
      "  'Variance Ratio Test Statistic (k=20) Heteroscedasticity Assumption': 0.0007646398954355974,\n",
      "  'Variance Ratio Test Statistic (k=20) Homoscedasticity Assumption': 0.7655801985572465},\n",
      " {'Variance Ratio (k=30)': 1.0054447472916037,\n",
      "  'Variance Ratio Test Statistic (k=30) Heteroscedasticity Assumption': 0.0008819254471013547,\n",
      "  'Variance Ratio Test Statistic (k=30) Homoscedasticity Assumption': 0.8829960534693014},\n",
      " {'Variance Ratio (k=40)': 1.007383025302277,\n",
      "  'Variance Ratio Test Statistic (k=40) Heteroscedasticity Assumption': 0.001029021845184675,\n",
      "  'Variance Ratio Test Statistic (k=40) Homoscedasticity Assumption': 1.0303005120741011},\n",
      " {'Variance Ratio (k=50)': 1.0086502431826903,\n",
      "  'Variance Ratio Test Statistic (k=50) Heteroscedasticity Assumption': 0.0010741842833486778,\n",
      "  'Variance Ratio Test Statistic (k=50) Homoscedasticity Assumption': 1.0755809312730416},\n",
      " {'Variance Ratio (k=100)': 1.0153961901671607,\n",
      "  'Variance Ratio Test Statistic (k=100) Heteroscedasticity Assumption': 0.0013415126178608341,\n",
      "  'Variance Ratio Test Statistic (k=100) Homoscedasticity Assumption': 1.3434284573260966},\n",
      " {'Variance Ratio (k=200)': 1.015704654116103,\n",
      "  'Variance Ratio Test Statistic (k=200) Heteroscedasticity Assumption': 0.0009639238446200716,\n",
      "  'Variance Ratio Test Statistic (k=200) Homoscedasticity Assumption': 0.9653299929053236},\n",
      " {'Variance Ratio (k=500)': 1.018216620766853,\n",
      "  'Variance Ratio Test Statistic (k=500) Heteroscedasticity Assumption': 0.0007055684744355339,\n",
      "  'Variance Ratio Test Statistic (k=500) Homoscedasticity Assumption': 0.7065863036900603},\n",
      " {'Variance Ratio (k=1000)': 1.0187822241562867,\n",
      "  'Variance Ratio Test Statistic (k=1000) Heteroscedasticity Assumption': 0.0005140701392295625,\n",
      "  'Variance Ratio Test Statistic (k=1000) Homoscedasticity Assumption': 0.5147582201029187}]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "# Generate random steps with mean=0 and standard deviation=1\n",
    "steps = np.random.normal(0, 1, size=1000000)\n",
    "# Set first element to 0 so that the first price will be the starting stock price\n",
    "steps[0] = 0\n",
    "# Simulate stock prices, P with a large starting price\n",
    "P = 10000 + np.cumsum(steps)\n",
    "# Test\n",
    "data = pd.DataFrame(P, columns=['Price'])\n",
    "result = estimate(data)\n",
    "pprint(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tbquant",
   "language": "python",
   "name": "tbquant"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
